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Abstract. We give an explicit description, via analytic subordination, of free 
multiplicative convolution of operator-valued distributions. In particular, the 
subordination function is obtained from an iteration process. This algorithm 
is easily numerically implement able. We present two concrete applications of 
our method: the product of two free operator-valued semicircular elements 
and the calculation of the distribution of dcd + (Pcd? for scalar- valued c and 
(i, which are free. Comparision between the solution obtained by our methods 
and simulations of random matrices shows excellent agreement. 



1. Introduction 

Free probability is a quite recent theory that has gained interest in the last 
few years. One of its main applications is on the field of random matrices. More 
specifically, it provides a conceptual way of understanding the distribution of the 
eigenvalues of several large random matrices. The variety of matrix models where 
free probability can be used is growing in accordance to the developments of the 
theory. 

The crucial requirement that allows the treatment of a random matrix model 
with the algebraic and analytical machinery of free probability is that the matri- 
ces involved should satisfy (asymptotically, as their size tends to infinity) freeness 
relations. Some of the most important random matrices, such as Wigner and Haar 
Unitary matrices, were shown, starting with the basic work [27 , to have these free- 
ness requirements among themselves and with respect to deterministic matrices. 

The applicability of free probability increased rapidly, in different directions, 
with the implementation of Voiculescu's operator- valued version of the theory. The 
main idea is that operator-valued freeness is a much less restrictive condition, but 
still most of the features of usual free probability theory are present. We are now 
able to work with block Gaussian matrices [TS", TT , rectangular random matrices of 
different sizes [5 and more complicated combinations of random and deterministic 
matrices |^, where scalar- valued freeness breaks up. 
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Whereas the analytic theory of scalar- valued free convolutions is far evolved (and 
thus we have now a variety of ways to deal with asymptotic eigenvalue distributions 
of matrices which are free), the same cannot be said about the status of operator- 
valued convolutions. In particular, at the moment we do not have an analytic 
description of operator- valued convolutions which would be easily and controllably 
implementable on a computer. Thus, numerical investigations for the above men- 
tioned random matrix models were done in a more or less ad hoc way. Our main 
aim is to improve on this situation and provide a coherent analytic description of 
operator-valued free convolution and show its usefulness for dealing with random 
matrix questions. In the present paper we concentrate on multiplicative free convo- 
lution, the question of additive free convolution will be addressed in another paper 

a. 

The present paper is motivated by the following problem (which was communi- 
cated to us by Aris Moustakas in this form in the context of wireless communications 
and, independently, by Raj Rao for the special case bi = 62 and ai = in the 
context of random graphs): If {ai,a2} and {61,62} are free, then it is not true in 
general that the elements ai^iai, a2b2a2 are free. This has made the distribution 
of c = aibiai + a2b2Ci2 quite inaccessible up to now. 

We observe, however, that the distribution of c is the same (modulo a Dirac 
mass at zero) as the distribution of the element 

/-jN f aibiai ^ a2b2a2 0\ _ /ai a2\ fbi 0\ /ai 0\ 

V 0/ vo vo W o; ' 

which in turn has the same moments as 

(2) f^' ^)='^AB. 

^ ^ \a2a1 J \0 62/ 

The advantage of this reformulation is that the matrices A and B are free with 
amalgamation over the algebra M2(C) of 2 x 2 constant matrices. Hence, the dis- 
tribution that we are looking for will be given by first calculating the M2(C)-valued 
free multiplicative convolution of A and B to obtain the M2(C)-valued distribution 
of AB and then getting from this the (scalar- valued) distribution of AB by taking 
the trace over M2(C). This has motivated us to look at the general problem of 
dealing with operator- valued free multiplicative convolutions. 

The standard way to deal with free multiplicative convolutions is through Voiculescu's 
S'-transform. The generalization of the S'-transform to the operator-valued situa- 
tion was found by Dykema [13 . (A certain version of the S'-transform, via a Fock 
space- type construction, appeared already in the work of Voiculescu [29].) Direct 
computations of operator- valued S-transforms for non-trivial elements are, how- 
ever, extremely hard. Approximations can be done on domains corresponding to 
domains of the Cauchy-Stieltjes transform which are located far away from the 
real axis. For practical purposes this is a problem, since we are interested in the 
behavior of the Cauchy transform close to the real axis, in order to recover the 
distribution by an Stieltjes inversion. 

In this paper we follow Biane's approach [9 to scalar-valued free multiplicative 
convolution via analytic subordination, which was later extended by Voiculescu 
j3m [3T| to the operator valued case. 
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Our main contribution to the theory is to find the subordination functions as 
iterative hmits, similar to what has been done in the scalar case [2 . We rely on 
the twisted multiplicative property of Dykema's ^-transform. 

This will allow us to set up fixed point equations to approximate effectively the 
values of the Cauchy transform of AB in the whole operator-valued upper half 
plane and in particular, close to the real axis. As inputs, we require the individual 
operator- valued Cauchy Transforms of A and B (or good approximations of these). 

We want to stress that in the operator- valued context there are rarely situations 
where one has an analytic formula for the involved Cauchy transforms. The best 
one can hope for are equations which determine those transforms. In order to be 
applicable to concrete problems one needs of course a way to solve these equations - 
in particular, to single out the correct solution; our operator-valued equations usu- 
ally have many solutions, only one of them corresponding to the wanted Cauchy 
transform. Since the characterization of the Cauchy transform among all solutions 
is by positivity requirements this is an intrinsic analytic characterization, which 
indicates that those equations cannot be solved by pure formal power series ex- 
pansion arguments. So what we need for a theory of operator- valued convolution 
which is also practically applicable is to set up our theory in such a form which 
can also numerically be implemented (and such that the theory provides arguments 
for the working of this implementation). What has become more and more appar- 
ent in the scalar- valued context, namely that the subordination formulation of free 
convolution seems to be the right choice, is even more prominent in the operator- 
valued context. Trying to solve an operator- valued free multiplicative convolution 
problem directly with the help of the operator-valued ^-transform becomes quite 
a challenging task very soon, whereas using the subordination formulation, as pre- 
sented in this paper, allows not only a satisfying analytic description of the involved 
transformations, but can also be implemented numerically very easily. 

In Section 2, we will present our analytic description, via subordination, of the 
free multiplicative operator-valued convolution. In Section 3, we will show the 
usefulness of our approach by implementing our method to calculate the distribution 
of the product of two operator- valued free random variables and by comparing this 
with histograms for corresponding random models. We consider there two different 
types of examples: (i) the product of two operator-valued semicircular elements, 
the individual Cauchy transforms obtained by the fixed point equations described 
in [14]; (ii) in the context of our original problem described above, we treat the case 
aba + a^ba^, where a is discrete and b is either discrete or a semicircular variable. 

2. Multiplication of operator- valued free random variables 

We will call an operator-valued non- commutative probability space a triple (A^, E, 5), 
where is a von Neumann algebra, B C is a VK*-subalgebra containing the 
unit of and E: Al — > 5 is a unit-preserving conditional expectation. Elements 
in A4 will be called operator-valued (or 5-valued) random variables. The distri- 
bution of a random variable x G A4 with respect to E is, by definition, the set of 
multilinear maps 

l^x := {^n- ^"""^ B: mn(6i,...,6n-i) =^xbixb2-'-xbn-ix],n G N}. 

We call the n*^ moment of x (or, equivalently, of /j^x)- It will be convenient 
to interprete mg as the constant equal to 1, the unit of B (or, equivalently, of A4) 
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and mf = E[x], the expectation of x. We denote by B{xi . . . the *-algebra 
generated by B and the elements Xi, . . . , in A^. 

Definition 2.1. Two algebras A\^A2 ^ M. containing B are called free with amal- 
gamation over B with respect to E (or just free over B) if 

E[xiX2 ■ --Xn] =0 

whenever n G N, Xj G Ai. satisfy E[xj] = and ij ^ ^j+i; 1 < j ^ ^ — 1- Two 
random variables x^y e M are called free over B if B{x) and B{y) are free over 
B. 

li x^y G M. are free over 5, then (ix+y and jixy depend only on ji^ and jiy. 
Following Voiculescu, we shall denote these dependencies by ji^ H l-^y and jix^ l^y^ 
and call them the free additive, respectively free multiplicative, convolution of the 
distributions fix and fiy. It is known [20l |29] that both ffl and M are associative, 
but M may fail to be commutative. 

2.1. Analytic transforms. A very powerful tool for the study of operator- valued 
distributions is the generalized Cauchy-Stieltjes transform and its fully matricial 
extension [29l[32: for a fixed x G M, we define G^(6) = E [(6 - x)-^] for ah 6 G 5 
so that 6 — X is invertible in M. One can easily verify that Gx is a holomorphic 
mapping on an open subset of B. Its fully matricial extension G^T^ is defined on 
the set of elements h G Mn{B) for which b — x In is invertible in M^(At), by 
the relation G^\b) = E (g) ldM^{c) [{b — x ^ ^n)~^] • It is a crucial observation of 
Voiculescu that the family {G'^^}n>i encodes the distribution fix of x. A succint 
description of how to identify the n*^ moment of x when {G'^^}n>i is known is 
given in [4 . 

In the following we will use the notation x > for the situation where x > and 
X is invertible; note that this is equivalent to the fact that there exists a real e > 
such that X > el. From the later it is clear that x > implies E[x] > (because 
our conditional expectations are automatically completely positive). 

From now on we shall restrict our attention to the case when x, y are self ad- 
joint, and (for most applications) nonnegative. In this case, one of appropriate 
domains for Gx - and the domain we will use most - is the operator upper half- 
plane EI+(5) := {b e B: ^b > 0}. Elements in this open set are all invertible, 
and W^{B) is invariant under conjugation by invertible elements in B. It has been 
noted in [3^ that G^T^ maps H+(Mn(5)) into the operatorial lower half-plane 
M-{Mn{B)) := -H+(Mn(5)) and has "good behaviour at infinity" in the sense 
that lim bG^J"^ (b) = lim G^^^ (b)b = 1. 

As, from an analytic function perspective, G^T^ have essentially the same be- 
haviour on H+(M^(5)) for any n G N, we shall restrict our analysis from now on 
to Gx '= G^x^ . However, all properties we deduce for this Gx-, and all the related 
functions we shall introduce, remain true, under the appropriate formulation, for 
ah n > 1. 

We shall use the following analytic mappings, all defined on H+(5); all trans- 
forms have a natural Schwarz-type analytic extension to the lower half-plane given 
by /(^*) = /(^)*; in all formulas below, x = is fixed in M.: 
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• the moment generating function: 

(3) = E [(1 - bx)-^ - 1] = E [(6-1 - - 1 = G,(6-i)6-i - 1; 

• The reciprocal Cauchy transform: 

(4) F,(6)=E[(6-x)-i]"'=G,(6)-i; 

• The eta transform (Boolean cumulant series): 

(5) Vx{b) = M/x(6) (1 + ^Ab)y' = 1 - 6E [(6-^ - x)-^] = 1 - 6F,(6-i); 

• In this paper we shall call this function "the h transform:" 

(6) h,{b) = h-\,(h) = - E [(h-^ - x)-'] = b-' - F,(h-^)- 

Based on the moment generating function, Dykema introduced the operator- 
valued version of Voiculescu's ^-transform [26 as an analytic mapping on Banach 
algebras (an earlier, less easily employed version can be found in [29]). It is easy to 
note that ^i,(6)(c) = E [(1 - hx)-^cx{l - 6x)-i] , so that ^^(0)(c) = E[cx] = cE[x]. 
Under the assumption that E[x] is invertible in B (and, in particular, when x > 0), 
the linear map ^^(0) becomes invertible, with inverse c ^ cE[x]-i, and so by the 
usual Banach-space inverse function theorem, has an inverse around zero, which 
we shall denote by - The S'-transform is defined as 

(7) 5^(6) = h-^il + 6)^-1(6), ||6|| small enough. 

Dykema showed [13, Theorem 1.1] that, whenever E[x] and E[?/] are both invertible 
in B, 

(8) S:,y{h) = Sy{b)S,{Sy{b)-^bSy{b)), \\b\\ SIRS^W CUOUgh. 

2.2. Three ways to the subordination function. Here we shall describe three 
ways to finding the analytic subordination functions for free multiplicative con- 
volution of operator- valued distributions. The analytic subordination has been 
proved in different contexts by Voiculescu and Biane [28i iHi i30| i31 . For the case 
5 = C, Biane showed that there exist analytic functions uji,uj2' C \ [0, +oo) 
C \ [0, +oo) which preserve half-planes and satisfy ^c{y) [{'^ ~ ^xy)~^] = (1 ~ 
uj2{z)y)~^, Ec{x) [{I - zxy)-^] = {1 - uji{z)x)-^. In particular, r]x{uji{z)) = 
r]y{uj2{z)) = r]xy{z)^ z e C \ [0, +oo). Voiculescu extended this relation, essen- 
tially in [30 , and made it more precise in [31], to the case of a general 5, under 
the assumption that A4 is endowed with a tracial state and E preserves this trace. 
Later, in [2 , Biane's subordination functions uji,uj2 were found as limits of an iter- 
ation process involving r]x and r]y. A precise description of such an iterative process 
for (very general) positive operator- valued random variables, in the spirit of [14], 
will be our main contribution in this section. 

Inspired by the shape of formula (|8|, we claim that for our purposes, the most 
appropriate writing of the operator-valued subordination phenomenon is the fol- 
lowing: 

Theorem 2.2. Let x > O^y = G A4 be two random variables with invert- 
ible expectations, free over B. There exists a Gateaux holomorphic map UJ2: {b e 
B: ^{bx) > 0} ^ IHI+(5), such that 

(1) riy{u;2{b))=rixy{b), ^{bx) > 0; 

(2) uj2{b) and b~^uj2{b) are analytic around zero; 



6 



S. BELINSCHI, R. SPEICHER, J. TREILHARD, AND C. VARGAS 



(3) For anybeB so that ^{bx) > 0, the map gi,: H+(5) H+(5), gb{w) = 
bhx{hy{w)b) is well-defined, analytic and 

uj2{b) = lim g^iw), 

n^oo 

for any fixed w G H+(5). 
Moreover, if one defines 0Ji{b) := hy{uj2{b))b, then 

Vxy{b) = uj2{b)r]^{uji{b))uj2^{b), ^{bx) > 0. 

Following the proof of this theorem, we shall mention several conditions under 
which some hypotheses, particularly the - rather inconvenient in practical applica- 
tions - invertibility requirement, can be weakened or entirely dropped. 

Proof We shall split our proof in several remarks, formulas and lemmas. For 
our purposes, a slight variation of the S'-transform will be more useful: we de- 
fine the sigma transform T^xi^) = b~^r]~^{b)^ again on a neighbourhood of zero. 
Elementary arithmetic manipulations show that T^x{b) = Sx{b{l — b)~^) and so 
^xy{b) = T^y{b)J^x{^y{b)~^bJ^y{b)). Using this relation we write on a neighbour- 
hood of zero in B: 

br^-^ib) = b'^^yib) 

= b^Ey(b)E,{Ey{b)-'bEyib)) 

= b'j:y{b){j:y{b)-'bEy{b))-'7J-\^y{b)-'bEy{b)) 
= bEy{b)7J-\Ey{b)-HEy{b)) 

= v;Hb)v^'{i%Hb))-'bv;Hb))- 

Now define uj2{b) = Vy^{Vxy{b)), again for ||6|| sufficiently small. We substitute in 
the previous relation rjxyib) for b to obtain 

r]y{uj2{b))b = r]xy{b)b 

= V;\Vxy{b))VxH{%HVxym-'Vxy{b)v;\v^^^^^ 

= ^2{b)rjx\^2{b)-'rjy{^2{b))u;2{b)). 
Recalling from equation (|6| the definition of the h transform, we obtain 

hy{uj2{b))b = ri~^{hy{ij02(b))uj2{b)), \\b\\ smah enough. 
An application of rjx on both sides gives r]x{hy{uj2{b))b) = hy{uj2{b))oj2{b)^ or 
(9) uj2{b) = bhx{hy{uj2{b))b) = gb{uj2{b)), \\b\\ smah enough. 

This shows us that UJ2 exists on a small enough neighbourhood of the origin, and 
is a fixed point for the map ^5 introduced in Theorem |2.2[ Moreover, this indicates 
that b~^uj2{b) is analytic around zero (a fact that follows quite easily also from 
the definition of UJ2 as 77"^ o rjxy). The most significant, however, is the fact that, 
under the conditions of analyticity of the two maps uJ2{b) and b~^uj2{b) around zero, 
equation ([9| uniquely determines a;2, and thus a function satisfying ([9| must also 
satisfy r]y o uj2 = rjxy^ and vice- versa. 

Now we observe that, under the additional assumption of the existence of a 
trace r on so that r = r o E, this function coincides with a function provided 
by Voiculescu: in |30j, Voiculescu proved that, whenever x and y are free over 
5, the range of the analytic map b ^ y ^ '^B{y) \(b — x — y)~^^ is included 
in B. {^B{y) denotes the conditional expectation with r onto the von Neumann 
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algebra generated by B and y.) It follows quite easily that the range of the map 
b ^ — ^B{y) [{y~^ — bx)~^] is also included in B. We claim that 

Since UJ2 up to this moment has only been defined on a neighbourhood of zero, 
we only need to verify the equality for ||6|| small. Indeed, the above relation is 
equivalent to ^B{y) [(1 — bxy)~^] = 1 — U02{b)y^ by the bimodule property of 
^B{y)- Inverting both sides of the equality and applying the conditional expectation 
E gives E [(1 — bxy)~^^ = E [(1 — uj2{b)y)~^^ . We invert, take a 6 as a factor and 
subtract 1 in order to get r]y{uj2{b)) = r]xy{b). As it was noted by Voiculescu 
[30 , and we shall argue below, y~^ — bx is invertible whenever ^{bx) > 0. Since 
the set {b e B: ^{bx) > 0} n {b e B : \\b\\ < e} is open for any e > 0, and 

lim ly~^ —KB{y) [{y~^ — ) = 0? it follows, by analytic continuation, 

||6||^0 V / 

that there exists an £ > so that 

UJ2: {beB: ^{bx) > 0} U {b e B : \\b\\ < e} ^ B 

(10) U2{b) = y-' - Esiy) [{y-' - bx)''] . 

Until now we have argued that an analytic map UJ2 satisfying parts (1) and (2) 
of our theorem exists and is unique (i) on a neighbourhood of zero for any B- 
valued non-commutative probability space (A1,E,5), and (ii) that this uj2 extends 
analytically to {b G B: ^{bx) > 0} when A4 is endowed with a tracial state r so 
that r o E remains a trace. We go now to the third way of identifying 002^ namely 
as a fixed point of an analytic mapping, method which will allow us to extend 
UJ2 to a set of the form {beB: ^{bx) > 0}U {b e B: \\b\\ < e} for any type of 
non-commutative probability space (A^,E, 5). 

Lemma 2.3. Assume that x > is invertible in M and b e {b e B: ^{bx) > 0}. 
IfceBisso that > 0, then ^h^ic) > E > 0. 

Proof. For simplicity, let us replace c G H+(5) by —c~^ G EI+(5), so that 

h,{-c-^)=E[{c^bx)-^]~^ -c. 

We shall split our problem in two and use the same method as in [T: assume that 
(p is an arbitrary positive linear functional on B so that (p{l) = 1. We define 

: C+ ^ C+, U{z) = p(e [{^{c + bx) + z^{c + bx))-^] . 

As > 0, and ^z are all strictly positive, it follows that ^{^{c-\-bx)-\-z^{c-\- 

bx))~^ < 0, and, since E is positive and faithful, ^E [(3?(c + bx) + z^{c + bx))~^] < 

0. Thus, ^ ^E [(3?(c + bx) + z^{c + > 0, and, in particular, invertible. 

Since (p > 0,(p{l) = 1, the Cauchy-Schwarz inequality tells us that '^f^p{z) > 
whenever > 0. We take 



lim = I lim E 



^{c^bx) 

h ^(c + bx) 



p(E[{^{c^bx))-^] 



8 



S. BELINSCHI, R. SPEICHER, J. TREILHARD, AND C. VARGAS 



The Nevanlinna representation [1, Chapter III] ahows us to write the inequahty 
^U{z) >(p(e [{^{c^bx))-^]~^^ for ah z G C+. Equahty holds at one point 

of C+ if and only if fcp{z) = cp (k [(^(c + bx))~^] z + real const ant. Taking z = i 
in the above we get 

^ip{E [{c^bx)-^]~^) > ip{E [{^{c^bx))-^]~^). 
Since this inequality holds for all states cp on we conclude that 

[(c + bx)-^] > E [(^(c + bx))-^] . 
We shall now prove that E [(^(c + bx))~^] ^ > ^c: 

E [(^(c + bx))-^] > ^ E [(^c + < i^c)-^ 

< 1 



^c{^c^^{bx)) 
1 



1 + (y^^ ^ ^{bx) (y^^ j < 1. 



< 1 



The last inequality is trivially true by functional calculus and the invertibility of 
^{bx). Putting the inequalities together gives 

ip{^E [{c^bx)-^y) = ^if{E [{c^bx)-^y) > ifi^c). 

Since this is true for all positive we get that h^x maps H+(5) into itself. 
Next, we make use again of the same trick: we define 

: C+ ^ C+, U{z) = ^(e [{c + ^{bx) + z^{bx))-^] - c) . 

As before, whenever > (the case c = c*, for example, is not excluded 
here), ^E [(c + -\- z^{bx))~^] — > for any z G C+. We take again 

limit as z ^ oo of f^p{z)/z to obtain ^E The same argu- 

ment used above implies that ^f^{z) > (p{E [{^{bx))~^] ^)^z, and so, for z = i 
we obtain ^ip (e [{c ^ ^{bx) ^ i^bx))-^y - > (p{E [{^{bx))-^y), for ah 
^ ^ O^V^ll) = 1- This implies ^h^xic) > E [(^(^x))""*^] as claimed. □ 

By the definition of the h transform ([6|, we have bhx{cb) = hbx{c). The previous 
lemma allows us to write 

^hy{w) > =^ ^hbx{hy{w)) > E > 0, 

i.e. ^gb{w) > E > for any w e H+(5), and b with ^{bx) > 0. 

Thus, we have shown that lands strictly into W^{B) whenever ^{bx) > and 
^hy{w) > 0. 

Remark 2.4. It has been shown in |4] that ^E [{w — y)~^] ^ > for any w G 
H+(B) andy = y* e M. By noting that -w~^ e H+(5) if and only if w G W^{B), 
we obtain that ^hy{w) > and ^hy{—w~^) > for any w G W^{B). 
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We can now improve on our previous statement: for x > 0, = and b e B so 
that ^{bx) > 0, 

> =^ ^hy{w) > ^hba:{hy{w)) > E [{^{bx))-^]~^ > 0, 
i.e. ^gb{w) > E [{^{bx))-^]~^ > for any w G IHI+(5), and b with ^{bx) > 0. 

Remark 2.5. The functions h associated to selfadjoints y G Ai have convergent 
power series expansions around the origin. Indeed, 

E[l -^^] -E [(1 -^^)"^]"^ 

= E [1 - wy] (e [(1 - wy)-'] -E[l- wy]''^ E [(1 - wy)''] 

= E[l -wy]w I f^E [y{wyr] - E[^]E[^^]- ) E [(1 - wy)''] . 



\n=l 



Thus, 
hy{w) 



(l-E[(l-^^)-i]-') 



= w 



E [1 — wy] w 



; E [y{wyr] - ny]E[wyr E [(1 - wy)-^] 



{1-E[yw]) 



E [yiwyT] - E[y]E[wyr ) E [(1 - wy)-'] ' 



\n=l 



■E[y] 



which gives the power series espansion of hy around zero and shows that hy{0) = 

We note that, as shown in the above Remark, there exists a £ > so that hy 
are defined on {w G B: \\w\\ < e} and 

max{||/.,H||, \\hy{w)\\ :weB, \\w\\ < e} < 2{\\x\\ + ||^||). 

Choosing be B with ||6|| < (e/41) • (||x|| + ||^||)"^ guarantees that \\hy{w)b\\ < e/10 



and \\bhx{hy{w)b)\\ < e/2, so g^ maps {w e B: \\w\\ < e} into {w e B: \\w\\ < 
s/2}. Thus, by the Earle-Hamilton Theorem [12, Theorem 11.1], gb has a unique 
attracting fixed point in {w G B: \\w\\ < e} which we shall denote by 0J2{b), and 

lim 56°"H ■■=uJ2{b) 

n^oo 

exists for all w e {w e B: \\w\\ < s}, \\b\\ < (e/41) • {\\x\\ + The corre- 

spondence b oj2{b) is clearly analytic, being a uniform limit of analytic maps. 
Moreover, on a small enough neighbourhood of zero, r]y{uj2{b)) = r]xy{b) by the 
uniqueness of the fixed point guaranteed by the Earle-Hamilton Theorem. 

For our fixed x^y G A4 given in the statement of our theorem, let us fix an £: > 
as above. Consider the set 



beB: 



< 



41 k 



ll^ll 



99 b- 



ll^ll 



This set is clearly open and connected (we can find the element b = i 82||a;-i||Q|a;|| + ||?y|| 
in both open sets whose union we considered). The above indicates that 
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(1) g: {b&B: \\b\\ < iTomw} >< > ^ B: \\w\\ <e}^{wGB: \\w\\ < §} 
given by gb{w) = hbx{hy{w)) is analytic; 

(2) g: [h e B: %hx) > 99||.-Miakll+||»y||) €B:^w> 200||.-i|K||.||+||,||) } 

{w G B: '^w > K[{^{bx))~^]~^} given by gb{w) = hi)x{hy{w)) is ana- 
lytic. 

The union of the two domains of g mentioned above is again connected, as we 
immediately note by identifying the point {b,w) = ( 82|^^^rfRHMD ' ^) ^^^^ 
sets. Since the following chain of implications holds, 

^bx) > n ^ {^{bx))-^ < - ^ E[{^bx))-^] < - ^ E[(^(6x))-^]-^ > 
the strict inclusion 



{w eB:^w> E[{^{bx))-^]-^} C 1^ G 5: 



200\\x-m\\x\\^\\y\\)} 

holds whenever ^(bx) > oon -mn ^ • 

Next, for any positive linear functional cp on B and b G B so that ^{bx) > 
99||a:-i||dkll + ||?y||) ^ ^^^^^ ^^^^ ^^^^ {^(^r(^))}nGN IS bouudcd. The argument 
uses the fact that the set {b e B: ^{bx) > 0} is convex, we choose bi G 5, 
^(^1^) > 99||x-i||(||.|| + ||^||) ^ ll^ill < 41(11.11 + 11^11) and consider t ^ ^{9tbHi-t)b,M)' 
The map [0,1] 3 t glb+(^i_i-^bS'^) lands, for a fixed w G EI+(5), entirely in 
H+(5) + ^ 99||a;-i||Q|a;|| + ||jy||) ? independently of n G N. Thus, there exists a small 
enough simply connected complex neighbourhood V of [0, 1], which does not depend 
on n so that V 3 t ^ (^) stih lands in EI+(5) for ah n G N. We obtain 

that all maps in the family 



[V3t^v{g°^^^,_,^,^{w))]^ 



take values in C+. This means that the family is normal (as a family of functions 
between complex domains). On the other hand, for |t| very small, we know that 

lim (w) = 0J2{th + (1 - t)W) e B, 

which, in particular, means that lim^^oo ^{glbj^(^i_t)hS^^^ ~ V^(^2(^^ + (1 ~ ^)^i)) 
exists and is finite. This, together with the above argued normality, implies that 
lim^^oo ('^)) exists as a holomorphic function from the given neigh- 

bourhood V of [0, 1] to C+. Now any linear functional on B has a unique Jordan 
decomposition as a linear combination of four positive linear functionals. Thus, for 
any if in the dual of 5, the family {\^{glll'j^(^i_^-^^_^{w))\ : t G n G N} is bounded. 
The uniform boundedness principle (see, for example, [23^ Lemma 1]) guarantees 
that {\\glJ^^(^-^_^^^^^{w)\\\ t G n G N} is bounded. The fact that 5 is a von Neu- 
mann algebra implies that {gl^{w)}nen niust have a w-convergent subsequence. 
Since ^g^'^{w) > E [(^(6x))~^] ^, it is clear that any such limit point must belong 
to H+(5). However, as noted before, lim^^oo ^(^^'^^(i-^)^^ ('^)) exists for all t eV^ 
^ B* . Thus, there can only be one limit point, i.e. lim^^oo 5'6^('^) = ^2(^) niust 
exist for all b with ^(6x) > 0. 

Note that in fact we have shown more: the limit function UJ2{^) •= linin^oo ^5^('^) 
is Gateaux holomorphic [TTl Definition 2.1] on all of the set {b ^ B: ^(bx) > 0}. It 
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is known to be holomorphic close to the origin of where r]y{uj2{b)) = rjxyih). Us- 
ing the same convexity trick as in the previous paragraph and the identity principle 
for the usual (scalar) analytic functions, we conclude that 

r]y{ijJ2{h)) = r]xy{h) and gij{u2{b)) = c^2(^), ^{bx) > 0. 

Thus, we have proved parts (l)-(3) of our theorem. The property of ui is trivial. □ 

The above theorem, as stated, has the inconvenience that it does not cover 
the posssible case of non-invertible positive x and the case of non-invertible IE[^]. 
However, when B is finite dimensional (a matrix algebra), we can use normality 
of some families of analytic maps from to itself to considerably improve our 
result. We shall use the same notations as in Theorem [2^ 



Proposition 2.6. Let B be finite- dimensional. For any x > 0, y = free over B, 
there exists a domain V C B containing and an analytic map 002'- V ^ H+(5) 
so that 

r]y{uj2{b)) = r]xy{b) and gb{uj2{b)) = uJ2{b), beV. 
Moreover, 0J2{b) = lim^^oo ^5^('^) for any w G W^{B), b gV. 

Proof We shall take as P = int{6 G B: ^{bx) > 0, 6 invertible in B}. By Remark 
2^ gb : EI+(5) EI+(5) is well-defined and analytic for any b with ^{bx) > 0. The 



existence of an attracting fixed point for ^5 when b is very close to the origin follows 



by exactly the same argument as in the proof of Theorem 2.2 Since the family {V 3 
b ^ 9b^{'^)}neN is normal, we conclude as before that 0J2{b) := lim^^oo ^5"^('^) 
exists and is analytic. (Here the fact that dim{B) < 00 is essential!) By the 
identity principle, it follows that U2{b) = gi){uj2{b)) for all 6 G P, as the relation is 
known to hold for b of small norm. 

Up to this point, we have shown the existence of an uj2 defined on the open set 
V D C+ • 1 which satisfies gb{uj2{b)) = 0J2{b)^ independently from the invertibility 
of either x or Since dim(5) < 00, the spectrum of E[?/] is a finite set in R, 

so it follows that non-invertibility of K[y] is equivalent to ]E>[y] G B having zero as 
eigenvalue. For any e > 0, E[?/ + e • 1] = E[y] -\- e ■ 1 is then invertible in B, and 
so is X + £ • 1 in M. If gh.ei"^) = h^x-\-s){hy-\-£{w)) ^ w G EI+(5), ^{b{x + e)) > 0, 
it follows that g^^^ ^5 uniformly on compact subsets of H+(5) as e ^ 0. Since 
for b of small norm, uj2{b) is an attracting fixed point for ^5, it follows that the 
small attracting fixed points of ^5^^, which we shall call 0J2{b)^ converge to uj2{b). 
Normality allows us to conclude that co'l ^ uj2- 

On the other hand, as seen in Theorem 12.21 



Since rjy and t]^^,^^)^^^^) rjxy^ we conclude that 

rjyOuj2 = r]xy on V. 

□ 

3. Numerical Implementation 



We now consider some of the practical implications of Theorem 2.2 in the com- 
putation of the free multiplicative convolution of operator- valued random variables 
- in particular, those which are represented by n x n matrices. The general frame for 
those examples will be the following. Our problems will be given (directly or after 
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some manipulations) by two matrices x = {xij) and y = (i/ij) where the entries of 
those matrices are hving in some non-commutative probabihty space r), and 
where, with respect to r, the entries of x are free from the entries of y. 

Let us set := Mn{A) = Mn{C)^A and consider the trace cp := tr(g)r : M ^ C 
and the conditional expectation E := id r : ^ M^(C). 

Thus X and y are elements in the non-commutative probability space {M,(p) 
and we are actually interested in the scalar- valued distribution of z := xy (or 
some variant thereof) with respect to = tr (g) r. However, the freeness between 
the entries of x and the entries of y with respect to r does in general not ensure 
freeness between x and y with respect to tr r. What it implies is operator- valued 
freeness between x and y in the non-commutative operator- valued probability space 
(A^, E, M^(C)). Thus we can first calculate the operator- valued distribution of z 
with respect to E = id r by using our free convolution results and then get the 
scalar-valued distribution of z with respect to = tr (g) r by applying the trace. 
More specifically, if Gz is the M^(C)-valued Cauchy transform of then tT{Gz) is 
the scalar-valued Cauchy-transform of z. Thus the spectral distribution of z with 
respect to tr (g) r is given, by virtue of the Cauchy-Stieltjes inversion formula, by: 

(11) d/j.{t)= lim ^^(tr(G^((t + ze)/n))), t^R 

where In is the nxn identity matrix and tr := -Tr is the normalized trace on n x n 
matrices. So our objective requires that we compute the Cauchy transform of xy 
at points of the form z/^, where z e C. 
Note now that by Q and (|6|, we have: 



(12) 



SO computing Gxy is equivalent to computing hxy. Moreover, by Theorem 2.2 and 
([5| we see that 



(13) Zhxy{zln) = UJ2{zIn)hy{uJ2{zIn)) 



Theorem 2.2 expresses the function uj2 as the limit of iteratively composing the 
function ^5 with itself. (Since in our case B = Mn{C) is finite-dimenisonal, we 
are actually in the realm of Prop. |2.6[ and thus do not have to bother about 
the invertibility assumptions.) It is now clear that given the operator- valued h 
transform of x and y, we can numerically compute the hxy, and thus the spectral 
distribution of the product xy. We consider several concrete examples below, in 
Sections 13.11 and 13.21 

One obstacle to applying this technique to general problems is the difficulty in 
analytically computing the h transform of many elements. It is easy to find the 



exact expression for the h transform of discrete distributions (see Section 3.2), and 
numerical methods exist to compute the h transform in other cases (see Section 3.1 
for example). 



3.1. The Product of Two Free Operator- Valued Semicirculars. Let 5i, 52, 

53, and 54 be free, semi-circular random variables, in some scalar- valued non- 
commutative probability space {A^ r) . Consider the matrices and 5*2 defined 
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by: 

(14) S, = h ''),S,= K+'' 

^ ^ \Sl S2J \ 2S4 S3 - 3S4 

Matrices and 5*2 represent limits of random matrices, where 5i,...,54 are 
replaced by independent Gaussian random matrices. 

As before, we set A4 := M2{A) = M2(C) and consider the trace cp := tr(g)r : 
^ C and the conditional expectation E := r id : ^ M2(C). 

We wish to compute the spectral distribution of (^'2 + cl2)Si in the scalar- 
valued probability space where c is some constant chosen large enough 
to make 5*2 + c/2 positive. Since Si and S2 are not free in {M^(p) we cannot 
invoke usual free probability theory to achieve our goal. However, Si and 5*2 are 
free operator- valued semicircular elements in the operator- valued probability space 
(A^, E, M2(C)). Thus we can do the calculations on the operator- valued level and 
in the end go down to the scalar- valued level by taking the trace. 

The first task is to compute the operator- valued h transforms, hs^ and hs^- In 
this case, an analytic equation for the h transforms is difficult to achieve. However, 
we can compute these h transforms numerically using the method described in 
[14]. In brief, this involves expressing the Cauchy transform of the operator- valued 
semicircular in terms of the fixed point of a contraction mapping. Specifically, if 
we define 

(15) W{h) = lim T^iWo) 

n^oo 

where Tb{W) = {-ib ^E[SbS]y^ , then Gs{b) = -iW{b). Note that we require 
the initial state Wq to satisfy ^(Wo) > 0; convergence of the above iteration scheme 
is ensured by arguments from [14 . In our case, with b = (^ij)f^j^i, we have 

E[SibSi]=(^''^l''^l''^^'' 

^ ^ V + ^11 + b22 

and 

'2611 + 2621 + 2612 + 4622 2611 + -2612 + 4621 - 6622 



HS2bS2, V^^ii + 46i2 - 2621 - 6622 46n - 6612 - 6621 + IO622 

We compare the spectral distribution of Si and ^2 computed using this method 
and the Cauchy- Stieltjes inversion formula to random matrix simulations in Fig. [l] 

Finally, using the numerically computed h transforms of Si and S2 + c/2 we used 
the iterative method discussed here to compute the h transform of their product. In 
Figure [1] we compare the distribution computed using our method to random ma- 
trix simulations of the ground truth spectral distribution of ^5*2+ cl2Si\/ S2 + c/2- 

For the sake of variety, we consider another operator- valued semi-circular exam- 
ple. Let now {si}^^^ be free semi-circular elements, and let: 

(16) 

(— lOSi 2^2 SOSsX / — 2^4 + 3S6 3^5+3056 Sq 

2s2 — 4s3 bsi and S2 = 3^5 + 30^6 54 + 55 + sq S4 

3O53 5Si IQSij \ Sq S4 40S4y 

We follow the same pattern as previously: applying the numerical method pro- 
posed in [14^ to compute the individual h transforms of S[ and (see Figure 
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Figure 1. Spectral distribution of Si (left), S2 (middle), and 
{S2 + 8. 5/2) 5*1 (right) - random matrix simulations (histogram) 
compared with numerically calculated density, using fixed point 
method of \JA^ for Si and 5*2 and using our method for (^2 + 
8.5/2)^1 . 



[2|, and then using our iterative method to compute the spectral distributions of 
(5^ + 85l3){S[ + 40/3) and {S!^ + 85l3){S[ + 75/3) (see Figure|3|. 




Figure 2. Spectral distribution of S[ (left) and ^2 (right) - ran- 
dom matrix simulations (histogram) compared with numerically 
calculated density using fixed point method of [14 . 




Figure 3. Spectral distribution of (5^ + 85/2)(5'l + 4O/3) (left) 
and {S2 + 85/2) (^1 + 75/3) (right) - random matrix simulations 
(histogram) compared with numerically calculated density using 
our method. 

3.2. Distribution of dcd-\- d?cd?. We consider now a special case of the problem 
from the Introduction, namely of finding the distribution of dcd + d^cd^ ^ where c 
and d are free from one another. As noted in the Introduction 



(17) 



'dcd + cPc(f 0\ _ fd d^\ fc 0\ f d 0\ 
0) ~ [0 ) [0 c) 0) 
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has the same distribution as 



(18) 



c 
c 



d^ 



Now, since c and d are free, the problem of calculating the distribution of dcd-\-d'^cd'^ 
has been transformed into one which can be solved numerically using the iterative 
method proposed in this paper. It is crucial to note that since c and d are free, we 

have that and ^ are free over the matrices M2(C). 

In order to apply our iterative method from Section 2, we need as input the 
operator- valued Cauchy transforms (or the h transforms) of the matrices x and 
where 



For instance, we have that 
(20) G,{b)=E[{b-x)-^] 
and a similar formula for G., 



E 



1 



^22 



-^12 
611 - c 



_{bii - c){b22 -c) - 612621 V ~^21 

Ty. (Note that the entries of the 2x2 matrix b — x 
commute and thus the usual formula for matrix inversion applies.) Thus we need 
to be able to calculate quantities like 

r [[{bii - c){b22 - c) - bi2b2i]~^ {b22 - c)] 

in order to calculate Gxib). 

For instance, if we assume here that c and d are both discretely distributed, such 
expressions can be readily written down in analytic forms. Consider the concrete ex- 
ample where c is uniformly distributed with discrete support {0.4, 0.7, 1, 1.3, 1.5, 1.7} 
and d is uniformly distributed with discrete support {0.5,1,1.5,2,2.5,3}. In this 
case, the distribution of dcd^d^cd^ computed using the iterative method proposed 
in this paper is shown and compared to histograms in Figure |4j 

We also compute the distribution of dcd^d^cd^ where d is uniformly distributed 
over support {0.4, 0.7, 1, 1.3, 1.5, 1.7} and c is a shifted semi-circular element (Figure 

i- 





Figure 4. Comparison of distribution of dcd + d^cd'^ where c 
is uniformly distributed over {0.4,0.7,1,1.3,1.5,1.7} and d is 
uniformly distributed over {0.5,1,1.5,2,2.5,3} (left) and where 
c is a shifted semi-circle and d is uniformly distributed over 
{0.4,0.7,1,1.3,1.5,1.7} (right). 
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